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1. INTRODUCTION 

Multivariate models are under-represented in the 
literature on spatial statistics. There is a basic rea¬ 
son for this: univariate models are sufficiently com¬ 
plicated to keep us busy. Genton and Kleiber have 
done a fabulous job compiling and investigating the 
available models, with a focus on the important class 
of models that they, with collaborators, introduced. 
This paper gives a solid state of the art and points 
out just how many holes there are in the theory and 
practice associated with these fields. This gives us 
licence to point out some other holes and to suggest 
some important directions for the future. 

2. THERE IS POWER IN A SPECTRUM 

If we were to quibble about one thing in Gen¬ 
ton and Kleiber’s paper, it would be that we dis¬ 
agree over the extent to which the class of multi¬ 
variate GRFs has been categorized. Note that this 
is different from explicitly constructing valid cross¬ 
covariance functions! To wit, if a multivariate GRF 
has a spectral representation, the spectral represen¬ 
tation given in Section 1.2 completely characterizes 
the class of stationary multivariate random fields 
that admit an absolutely continuous spectral mea¬ 
sure. This represents a large chunk of interesting 
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GRFs. We note that the paper, by restricting the 
cross-spectral densities to be real, implicitly assumes 
that Cij{ h) = Cij(— h), when the minimal necessary 
requirement is only that C t j(h) = Cji{— h), which al¬ 
lows for phase differences between the model compo¬ 
nents. The representation can then be employed con¬ 
structively as follows. Let u —> S(u>) be a mapping 
from to the set of Hermitian nonnegative definite 
matrices, the elements of which the cross-spectral 
densities, denoted f t j in the paper, are here subject 
to /y(u>) = fji(uj). Then, for any complex, matrix¬ 
valued function L(u;) such that Lij( cj) = Ljj(— uj) 
and S(-) = L(-)L( : ), 

x(s) — / L(w)e <8 - w dW(w) 

(i) JkJ 

= fj 

J«. d JR d 

where cfW(-) € C p and dW(-) £ are Gaussian 
white noise processes on understood as ran¬ 
dom measures with dWi(uj) = dWi{—u), E[dW(w) • 

dW(w')] = 5{u - u/)I du, and E[dW(s) dW(s')] = 
5(s — s')Ids (Adler and Taylor (2007); Lindgren 
(2012)). This representation only covers multivari¬ 
ate GRFs with absolutely continuous spectral mea¬ 
sures; however, the same procedure applies to fields 
with an atomic spectral representation. The ab¬ 
stract feature that is hiding in all of this specificity is 
that we are explicitly constructing a square root of 
the multivariate covariance operator and using this 
square root to filter the multivariate white noise. On 
a compact domain, the covariance operator is a com¬ 
pact, trace class operator, and so this square root is 
well defined using the usual functional calculus. 

Another reason to further emphasize this spec¬ 
tral representation is that it is not only constructive 
in its own right, but also useful when transformed 
back to the nonspectral domain. Kernel convolu¬ 
tion methods (Higdon (1998)) have a storied his¬ 
tory in univariate spatial statistics and their gen¬ 
eralization to the multivariate case is straightfor¬ 
ward (Simpson, Lindgren and Rue (2012); Bolin and 
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Lindgren (2013)). Their advantage is that it is never 
necessary to identify the spectrum of the process 
or, in fact, the cross-covariance structure. Rather, 
for any L 2 matrix-valued function K(-, ■)i x ( s ) = 
j K(s, s') dW (s') is a valid mutlivariate GRF, which 
can be approximated by (carefully) approximating 
the corresponding integral with a sum. In fact, there 
is nothing special about white noise in this situation, 
W(-) can be any independently scattered R^-valued 
random measure. This leads to a natural way to 
construct non-Gaussian random fields (Aberg and 
Podgorski (2011)). 

Different choices of K(-, •) generating the same co- 
variance function will in the non-Gaussian case af¬ 
fect the dependence structure. Similarly, the choice 
of square root in the spectral representation be¬ 
comes relevant, and the two integrals in (1) are no 
longer guaranteed to give the same process model. 
We also note that K(s, s') does not have to be a func¬ 
tion of s — s' and, hence, this resulting field does not 
need to be stationary. 

An important advantage to the multivariate spec¬ 
tral and convolution kernel constructions is that 
the spectral operator L(u;) and the kernel matrix 
K (•, ■) can reflect the modeler’s knowledge about the 
physical process under consideration. This can lead 
to useful, informative covariance structures that 
are tailored to the specific application. This idea 
falls under the auspices of the physics-constrained 
cross-covariance specifications mentioned in Sec¬ 
tion 7.3, except for one key difference: while a cross¬ 
covariance structure is present, identifying it is not 
necessary for inference, either conceptually or com¬ 
putationally ! 

3. THINK LOCAL, ACT LOCAL 

In fact, there are many parts of the above con¬ 
struction that we can live without. The SPDE ap¬ 
proach introduced by Lindgren, Rue and Lindstrom 
(2011) as a computationally efficient reformulation 
of GRFs is an example of the same procedure where 
we never construct the kernel matrix. Instead, mul¬ 
tivariate models can be constructed using systems 
of equations, such that the kernel corresponds to 
(matrix-valued) Green’s function of some linear par¬ 
tial differential operator (Hu et al., 2013a, 2013b, 
2013c). Using the partial differential operator con¬ 
struction has several advantages over the direct 
kernel specification. First, when the representation 
is Markovian, it allows us to localize the process. 


This is especially convenient when moving to non- 
Euclidean spaces; the great tragedy of spatial statis¬ 
tics is that the Earth turned out not to be flat. The 
second major advantage is that this localization al¬ 
lows us to construct local approximations to the 
resulting random fields. For an n-dimensional ap¬ 
proximation to a field observed at N points, this re¬ 
duces the computational cost of fitting the field from 
0((pN) 3 ) to 0(pN + p 3 n 3 / 2 ), for the case d = 2. 
This also compares well with non-Markovian meth¬ 
ods, such as the aforementioned n-dimensional con¬ 
volution kernel methods that have computational 
cost of 0(pn 2 N + p 3 n 3 ). Given that we are now liv¬ 
ing in the age of “big data,” this is a serious ad¬ 
vantage to the SPDE specification. The third major 
advantage is that it is straightforward to construct 
nonstationary models by locally varying the partial 
differential operator in the model. This corresponds 
to the “physics” view, where dependency structures 
are specified locally and extended to a global covari¬ 
ance structure through a conditioning argument. In 
our experience, this is an extremely powerful tool 
for specifying useful covariance and cross-covariance 
structures. 

4. WITH LOW POWER COMES GREAT 
RESPONSIBILITY 

One of the principal challenges that we have en¬ 
countered when applying likelihood methods to mul¬ 
tivariate GRFs is that their likelihood surfaces tend 
to be flat. This is perhaps not a surprise. If fitting 
p univariate models requires the estimation of 0(p) 
parameters, then a p-conrponent multivariate model 
will require the same data to be informative about 
0{p 2 ) parameters. This parameter inflation becomes 
noticeable already for bivariate models, but for large 
p it is a serious issue, that can sometimes be par¬ 
tially alleviated by using a low rank model, ideally 
motivated by problem-specific knowledge. Another 
option is to impose a sparse structure on the lin¬ 
ear filter matrix operator or to impose constraints 
between the parameters. 

While flat likelihoods resulting from an exploding 
parameter space are annoying, there is a far more 
pathological problem for univariate GRFs. Funda¬ 
mentally, the range and variance parameters are 
not identifiable under the usual infill regime (Zhang 
(2004)). This leads to a ridge in the parameter space 
that can only be resolved using careful prior mod¬ 
eling (Simpson et al. (2014)). Ridges will also se¬ 
riously challenge numerical optimizers and MCMC 
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schemes unless they use enough second-order infor¬ 
mation to resolve it. It is currently unclear to what 
extent these problems extend to multivariate mod¬ 
els; however, we suspect that they do, due to the 
aforementioned parameter inflation. 
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